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Computational  Deformation  and  Shearbands  We  begin  by  summariz¬ 
ing  our  work  on  shearband  formation  [5];  the  second  sub-section  outlines  our 
work  on  adaptative  implicit-explicit  methods  for  elastic-plastic  deformation 

PI- 

Shearbands  The  equations  governing  the  motion  of  a  granular  continuum  are 
the  balance  laws  for  mass  and  momentum.  A  particularly  simple  model  to 
study  is  the  anti-plane  shearing  of  a  continuum.  If  (x,  y )  denote  the  in-plane 
coordinates,  z  the  direction  orthogonal  to  this  plane,  and  v  the  velocity  in 
the  z-direction,  then  the  continuity  equation  is  identically  satisfied  with  a 
constant  density,  and  the  only  non-trivial  momentum  equation  is  for  the 
z-direction, 

dtv  +  ce2dxTxz  +  ce2dyTyz  =  0.  (1-1°) 

T  is  the  Cauchy  stress  tensor,  where  compressive  stresses  are  taken  to  be 
positive. 

We  now  must  augment  (1.1a)  with  suitable  constitutive  information  de¬ 
tailing  the  evolution  of  the  stresses  Txz  and  Tyz.  To  this  end,  we  follow  Scha¬ 
effer  [10]  and  adopt  a  von-Mises  type  Yield  condition  with  strain-hardening, 
and  a  possibly  non-associative  flow  rule. 


a(<7  +  [J-  tf(<7,*)£<7]  Vu  =  0.  (1.16) 

Here,  a  —  (TXJ,  Tyz')'  and  the  *  denotes  transpose.  The  rotation  3?<r  gives 
the  direction  of  the  strain-rate  tensor,  where 

«  _  (  cos  (a)  sin(a) 
y  —  sin(a)  cos(a) 

and  a  is  the  angle  of  non-associativity.  Finally,  H(a ,  a)  contains  the  consti¬ 
tutive  information  related  to  loading  and  hardening: 


h(d)  +  |<x|2  cos(a) 

where  ct(x)  =  max0< 3<t  |<7(x,  s)j  is  the  maximum  previous  stress  level  (and 
is  an  equivalent  measure  as  the  total  accumulated  shear  strain  7).  A  typical 
model  for  the  hardening  function  is  h(a)  =  hoy/l  —  \<r\.  The  elastic  compo¬ 
nents  are  modeled  by  linear  elasticity.  The  switch  <  l  >  is  1  if  the  material 


is  loading  plastically  (i.e.  If  \cr(x,  f)|  =  |<7(x)|),  and  0  otherwise. 
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In  writing  the  equations  above,  we  have  non-dimensionalized  relative  to 
parameters  typical  in  a  biaxial  laboratory  test.  In  this  setting,  the  maximum 
possible  total  stress  \cr\max  =  1  and  the  elastic  wave  speed  ce  «  106.  In  order 
to  make  the  calculations  feasible  with  a  reasonable  expenditure  of  computing 
resources,  we  compute  with  an  artificial  wavespeed  of  1. 

It  is  convenient  to  examine  the  one-dimensional  system  to  get  a  feeling 
for  the  fundamental  dynamics  of  the  system.  If  we  assume  dy  —  0,  the 
resulting  system  is  hyperbolic  whenever  |oj  <  d„it ,  where  the  critical  stress 
&crit  depends  on  the  degree  of  non-associativity.  In  particular,  a„it  =  1  when 
a  =  0  (i.e.,  the  system  is  always  hyperbolic  when  the  flow  rule  is  associative) 
and  0^  <  1  when  a  ^  0.  Similar  statements  hold  in  the  two  dimensional 
setting. 

Roughly  speaking,  a  shearband  forms  whenever  any  one  component  of 
stress  reaches  a  maximum  and  begins  to  decrease,  even  if  the  total  stress 
(e.g.  the  squared  norm  of  the  stress  tensor)  is  still  increasing. 

Before  equations  (1.1)  lose  hyperbolicity,  the  system  consists  of  a  single 
conservation  law  for  momentum  balance,  and  two  stress-rate  equations  which 
are  not  in  conservation  form.  The  paper  [5]  details  a  Godunov-type  scheme 
which  integrates  the  governing  equations  numerically,  allowing  for  loss  of 
hyperbolicity  and  the  formation  of  shearbands.  The  Godunov  scheme  is 
related  to  ideas  in  (1,  11] 

Our  basic  scheme  consists  of  four  parts:  (i)  monotone  slope  determina¬ 
tion;  (it)  characteristic  tracing;  (tit)  wave  interaction  (tv)  conservative  up¬ 
date  of  momentum  and  time-centered  stress  update.  This  basic  scheme  must 
be  amended  wherever  hyperbolicity  is  lost.  In  that  case,  a  shearband  is  as¬ 
sumed  to  form  in  a  cell  which  experiences  change-of-type.  This  shearband  is 
treated  as  an  internal  boundary;  plausible  boundary  conditions  proposed  by 
Schaeffer  [10]  are  imposed  at  the  boundary  and  characteristic  tracing  is  used 
to  update  the  dependent  variables  on  each  side  of  this  boundary.  Integration 
in  the  rest  of  the  domain  is  not  affected. 

Results  of  our  computations  illustrate  how  the  loss  of  hyperbolicity  at 
one  location  promotes  unloading  nearby,  thus  localizing  deformation  into  a 
shearband.  This  unloading  provides  a  nonlinear  ‘regularization’,  in  that  loss 
of  hyperbolicity  does  not  lead  to  global  blow-up  of  the  solution.  Rather,  a 
smooth  solution  may  lose  differentiability  once  the  system  loses  hyperbolicity, 
but  the  strains  remain  piecewise  smooth  with  jump  discontinuities.  Further 
computational  experiments  are  in  progress,  allowing  for  closer  examination 
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of  the  phenomenon  of  unloading. 

Implicit- Explicit  Computations  As  mentioned  above,  the  speed  of  elastic 
waves  in  a  granular  material  are  on  the  order  of  106  cm. /sec.  In  compu¬ 
tations,  this  large  wavespeed  imposes  a  stringent  constraint  on  the  allowable 
timestep  for  an  explicit  calculation,  a  timestep  too  small  for  a  fully  explicit 
numerical  simulation  of  an  experimental  biaxial  test.  To  overcome  this  sta¬ 
bility  constraint,  we  are  developing  an  Adaptatively  Implicit-Explicit  (AIE) 
Godunov  method  for  wave  propagation  [7].  The  motivating  feature  is  this 
work  is  that  most  elastic  waves  are  weak  and  do  not  contribute  substantially 
to  the  dynamics.  The  most  important  waves  are  (1)  waves  of  strong  plastic 
response,  and  (2)  unloading  waves  from  the  shearband. 

Previously,  Fryxell  et.  al.  [4]  developed  an  implicit  Godunov  method  for 
gas  j  namics;  this  method  is  fully  second-order  in  time,  using  a  trapezoidal 
rule  integration  in  time.  The  price  paid  for  second-order  accuracy  is  that  the 
linear  algebra  system  to  be  solved  at  each  gridpoint  consists  of  2N  equations, 
where  N  is  the  number  of  conservation  equations  to  be  solved.  Further, 
the  desire  for  monotonicity  limits  the  timestep  to  approximately  twice  the 
explicit  timestep.  This  timestep  limit  is  still  too  small  for  our  elastic-plastic 
problem,  and  the  role  of  plasticity  in  our  system  requires  new  ideas. 

Following  ideas  of  P.  Colella  and  J.  Trangenstein,  we  have  begun  a  pro¬ 
gram  to  develop  the  AIE  method  for  elastic-plastic  deformation.  The  idea 
is  to  switch  smoothly  between  explicit  and  implicit  time-iDtegration  for  each 
wave,  independent  of  other  waves.  For  waves  which  are  computed  explic¬ 
itly,  the  method  proceeds  as  above;  for  waves  which  are  computed  implicitly, 
we  reduce  to  first-order  accuracy  in  time  (currently  using  backwards  Euler) 
and  space  (no  characteristic  tracing).  We  use  the  Godunov  methodology  of 
averaging  (in  space)  the  solution  of  Riemann  problems  to  ensure  good  disper¬ 
sion  characteristics.  While  the  backward  Euler  time  integration  may  dampen 
some  waves  which  should  not  decay,  this  integrator  does  give  a  monotone  (in 
time)  solution.  Further,  the  linear  algebra  problem  consists  of  N  equations 
at  each  gridpoint. 

To  date,  we  have  written  an  AIE  code  for  solving  the  elastic-plastic  model 
of  Antman  [1].  This  model  consists  of  2  conservation  laws  (in  1  space  dimen¬ 
sion),  plus  a  single  ODE  for  the  accumulated  shear.  The  fundamental  sim¬ 
plification  of  this  model  is  that  the  system  is  always  strictly  hyperbolic.  We 
have  made  comparisons  with  Trangenstein’s  analytic  solution  to  the  Riemann 
problem  for  this  model,  and  with  his  explicit  Godunov  calculations.  While 
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our  current  version  works  well  for  relatively  weak-wave  solutions  and  small 
multiples  of  the  explicit  timestep  (order  5),  performance  for  strong-wave  so¬ 
lutions  and  Vu)  large  timesteps  is  less  satisfactory.  We  are  currently  working 
on  an  iteration  scheme  to  overcome  these  limitations.  (Note:  Trangenstein’s 
explicit  computations  also  degraded  for  strong-wave  solutions,  and  he  mod¬ 
ified  the  wave  tracing  part  of  the  algorithm  to  obtain  better  performance. 
We  are  examining  how  his  modifications  may  be  adapted  into  our  method.) 

Our  primary  goal  in  developing  AIE  is  to  to  solve  the  anti-plane  shear 
model.  We  are  studying  the  modifications  necessary  to  adapt  our  current 
code  for  this  problem.  The  most  vexing  difficulty  is  how  to  accurately  (i.e. 
explicitly)  capture  unloading  waves  near  the  shearband  while  treating  other 
elastic  waves  implicitly.  Our  current  strategy  is  to  use  AIE  for  the  initial 
loading  of  the  sample,  until  the  maximum  stress  .is  large  (say  |<r|  ss  O-TS^crit), 
and  then  switch  to  fully  a  explicit  integration,  to  capture  shearbands  and 
unloading  waves  accurately. 

!  Gradient  Theory  Significant  work  is  required  to  understand  the  true  nature 
I  of  the  shear  band  before  physically  “correct”  jump  conditions  can  be  applied. 
Schaeffer’s  suggestion  in  [10]  is  to  apply  the  same  constitutive  relations  to 
the  rapidly  shearing  material  inside  the  band  but  include  a  new  (phenomeno- 
j  logical)  lengthscale  i.e.  the  band  thickness.  An  alternative  approach  is  to 
identify  the  the  most  important  dissipative  processes,  which  provide  a  high 
frequency  cutoff  to  the  strain-rate  blow-up.  Then  it  may  be  possible  to  de¬ 
rive  jump  conditions  in  the  limit  as  the  strength  of  these  dissipative  processes 
approaches  zero. 

We  have  studied  continuum  theories  which  provide  a  dissipation  mecha¬ 
nism,  with  the  intention  of  (i)  examining  well-posedness  of  these  models,  and 
(it)  deriving  jump  conditions.  To  this  end,  we  began  with  a  Cosserat  model  as 
a  possible  regularization  method.  The  couple-stresses  in  the  Cosserat  con¬ 
tinuum  do  provide  a  high-frequency  cutoff  and  regularization  for  shearing 
motions.  But  the  model  fails  to  provide  regularization  for  dilation. 

Therefore  we  have  examined  higher-order  gradient  theories  as  a  potential 
model  (see  [12]).  This  theory  consists  of  mass  and  momentum  balance,  plus 
constitutive  equations  which  can  be  written  as 

\devT\  -  p(7,A7)  (2a) 
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irV«  =  /9(7,A7)  .  (2c) 

These  equations  express  the  yield  condition,  flow  rule,  and  dilatancy,  depend¬ 
ing  on  the  accumulated  shear,  7,  and  the  Laplacian  A7.  The  accumulated 
shear  evolves  according  to  dt~l  =  \devV^\.  V  is  the  strain  rate, 

Vij  =  +  fyi) 

and  yM  denotes  the  plastic  part  of  the  strain-rate.  Also,  trV  denotes  the 
trace  of  V,  trV  =  Vu,  and  the  deviator  devV  =  V  —  (trV)I  where  /  is 
the  identity  matrix.  The  resulting  system  of  equations  differs  from  a  more 
conventional  Flow  Theory  of  plasticity  in  that  a  PDE  for  7  must  be  solved. 

The  analysis  of  simple  shear  is  completed,  and  analysis  of  general  flow 
is  partially  completed  [8].  In  the  case  of  simple  shear  with  a  velocity  u 
depending  on  the  single  space  variable  y,  the  governing  equations  may  be 
written 

1 

fa  ~  dy  fi{ 7,  dvyl)  =  0 

da  -  ^  dvu  =  0 

Linearizing,  looking  for  exponential  solutions  e’v*+*^t  and  solving  the  re¬ 
sulting  eigenvalue  problem,  we  find  A  =  —  fti,  where  /i  1  and  fx2 

are  the  partial  derivitives  of  fi  with  respect  to  its  first  and  second  argu¬ 
ments.  In  classical  Flow  Theory  of  plasticity,  is  positive  until  a  critical 
value  of  the  accumulated  shear  is  reached  and  then  becomes  negative;  this 
non-monotonicity  is  the  source  of  ill-po?edness  and  shearband  formation.  If 
H2  <  0,  the  incorporation  of  the  higher  order  gradient  regularizes  the  equa¬ 
tions  and  leads  to  a  well-posed  system.  However  the  system  may  be  unstable 
for  a  range  of  wavenumbers  £,  and  concentrate  deformation  in  a  narrow  re¬ 
gion. 

Further  analysis  of  the  simple  shear  problem  also  shows  that  the  system 
exhibits  shearband-like  solutions,  similar  to  those  found  in  [2]. 


6 


\ 

\ 


\ 


After  completing  the  calculation  showing  well-posedness,  we  plan  to  ex¬ 
amine  whether  this  theory  can  be  used  in  an  anti-plane  shearing  problem,  as 
a  means  of  deriving  jump  conditions. 

Related  Systems  We  have  studied  two  issues  related  to  the  well-posedness 
of  elastic-plastic  deformation  problems.  The  first  issue  concerns  constitutive 
relations  with  rotational  symmetry;  the  second,  a  relaxation  phenomenon. 
Rotational  Symmetry  The  issue  of  isotropy  within  hyperelastic-plastic  models 
in  multiple  dimensions  is  important.  In  particular ,  if  the  stresses  are  assumed 
to  be  iotationally  invariant,  then  the  governing  equations  are  not  strictly 
hyperbolic.  The  origin  in  state  space  is  an  umbilic  point.  Existence  of  an 
umbilic  point  becomes  a  serious  issue  when  motions  which  include  the  origin 
are  permitted.  In  particular,  the  usual  high  order  techniques  for  solving 
hyperbolic  conservation  laws  fail  for  such  motions.  With  H.  Freistuhler  we 
are  examining  this  issue  [3]. 

For  rotationally  invariant  materials,  the  stress  is  usually  taken  to  be  a 
function  of  the  strain  of  the  form 

t  m  =  m\2)v- 

We  examined  a  model  of  this  kind,  numerically  solving  the  Riemann  problem 
for  the  2x2  system 

dtU  +  dx{\U\2U)  =  0  (3) 

where  U  =  (u,  u).  (Remark:  There  is  a  precise  isomorphism  between  the 
wave  patterns  of  (3)  and  (a  specific)  part  of  the  wave  pattern  of  any  generic 
rotationally  degenerate  system  of  hyperbolic  conservation  laws,  including 
isotropic,  neo-Hookian  elasticity.)  It  is  unclear  what  additional  effects  plastic 
yielding  will  have  on  the  computational  solution. 

Equation  (3)  is  hyperbolic,  with  eigenvalues  Aa  =  |17|2  and  Ar  =  3|f7|2, 
and  right  eigenvectors  ra  =  (-u,u)  and  rT  =  (u,v).  Note:  the  slow 

azimuthal  field  is  linearly  degenerate  in  the  sense  of  Lax,  while  the  fast  radial 
field  is  nonlinear  away  from  the  origin. 

We  solved  the  Riemann  problem  for  (3)  using  two  different  types  of  nu¬ 
merical  method,  viz.  a  random  choice  scheme,  and  a  Godunov  scheme.  The 
random  choice  method  is  only  first  order  accurate  in  space  and  time,  but 
introduces  no  artificial  dissipation  into  the  calculation.  In  contrast,  the  Go¬ 
dunov  scheme  we  implemented  is  second  order  accurate,  but  incorporates 
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a  selective  amount  of  diffusion  into  the  different  wave  families.  This  diffu¬ 
sion  corrupts  the  solution  for  special  Riemann  data.  Let  us  explain  this  last 
statement. 

For  Riemann  data  Ui  and  UT  which  is  colinear  and  lie  on  opposite  sides  of 
the  origin,  there  are  two  possible  constructions  of  a  solution.  For  a  specific 
example,  let  us  consider 

Ul  =  ( Ut,Vl )  =  (1,0)  Ur  =  (Ur,Ur)  =  (-1,0) 

One  may  construct  a  solution  which  follows  the  azimuthal  field  around  from 
(1,0)  to  (-1,0),  establishing  a  contact  discontinuity.  Also,  one  could  consider 
the  reduced  system 

dtu  v  dxu3  v  =  0 

The  solution  of  this  reduced  system  follows  the  construction  of  Wendroff  and 
Liu,  and  u  undergoes  a  composite  shock-rarefaction. 

These  different  constructs  point  to  a  non-uniqueness  in  the  solution  oper¬ 
ator  for  (3).  In  particular,  we  have  demonstrated  in  [3]  that,  for  nearly  colin¬ 
ear  Riemann  data,  the  Godunov  scheme  approximately  mimics  the  reduced- 
system  construction,  while  the  random  choice  scheme  (by  design)  mimics  the 
centered  wave  construction. 

We  believe  this  result  illustrated  the  difficulties  to  be  overcome  in  the 
numerical  simulation  of  hyperelastic  material  models.  This  project  is  contin¬ 
uing,  as  Freistuhler  and  I  study  other  numerical  algorithms  for  rotationally 
invariant  systems,  in  particular,  questions  of  convergence  for  the  Cauchy 
problem. 

Viscoelastic  Relaxation  Work  is  in  progress  studying  relaxation  mechanisms. 
Consider  the  system 


dtu  —  dxv  =  0  (4a) 

dtv-dxP  =  0  .  (46) 

If  we  provide  constitutive  information  P  =  P(u,v),  the  system  is  closed  and 
is  used  as  a  model  of  deformation  of  an  elastic  bar.  With  P  =  Pre/  =  u3/3— u, 
the  system  is  hyperbolic  whenever  |u|  >  1  and  elliptic  otherwise.  This  model 
is  used  to  study  change-of-phase  in  elastic  materials.  To  prove  existence 
of  solutions,  Eq.  4a-4b  must  be  augmented  with  an  admissibility  criterion 
describing  the  types  of  discontinuities  allowed  in  the  solution.  A  result  of 
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Shearer  shows  non-uniqueness  of  the  solution  for  certain  initial  data.  That 
is,  different  solutions  emerge  for  initial  data  in  a  special  region  in  the  ( u,v )- 
plane.  One  solution  consists  of  two  rarefaction  waves  coupled  with  two  phase 
jumps,  and  is  obtained  in  the  limit  when  a  viscosity  or  viscosity-capillarity 
regularization  is  used. 

We  are  examining  a  different  kind  of  regularization.  Assume  P  is  de¬ 
scribed  by  its  own  evolution  equation,  given  by  a  viscoelastic  relaxation  con¬ 
stitutive  model,  and  augment  Eq.  4a  and  4b  with 

6tP  +  E2dtu  =  ^  (4c) 

T 

where  r  is  a  relaxation  time.  In  what  sense  do  solutions  of  the  relaxation 
system  Eq.  4a-c  converge,  as  r  — »  0,  to  solutions  of  the  p-system  Eq.  4a-4b, 
with  P  =  PrS  As  part  of  his  Ph  D.  thesis,  Mr.  Yigong  Ni  is  studying  this 
problem  numerically  [9],  and  he  has  partial  analytic  results  on  convergence 
as  r  — ♦  0.  One  particular  numerical  finding  is  that  the  relaxation  limit 
solution  is  not  the  same  as  the  viscous  regularization  solution.  The  analysis 
gives  a  partial  answer  to  the  question  of  convergence  for  smooth  solutions, 
and  also  provides  asymptotic  behavior  near  discontinuities.  These  result 
will  be  important  in  our  granular  flow  work,  when  considering  viscoplastic 
constitutive  relations  (i.e.  relaxation  systems  such  as  Eq.  4,  with  plasticity 
included). 
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